Structure and Dynamics of Liquid Iron under 
Earth's Core Conditions 



On 



(N 



> 
On 

cn 
in 
o 
o\ 

c3 



i 

c 

o 
o 



13 



D. Alfe 1 , G. Kresse 2 and M. J. Gillan 3 

1 Research School of Geological and Geophysical Sciences 
Birkbeck College and University College London 
Cower Street, London WC1E 6BT, UK 
2 Institut fur Materialphysik, Universitdt Wien, Strudlhofgasse 4, A-1090 Wien, Austria 
3 Physics and Astronomy Department, University College London 
Cower Street, London WC1E 6BT, UK 

First-principles molecular dynamics simulations based on density-functional theory and the pro- 
jector augmented wave (PAW) technique have been used to study the structural and dynamical 
properties of liquid iron under Earth's core conditions. As evidence for the accuracy of the tech- 
niques, we present PAW results for a range of solid-state properties of low- and high-pressure iron, 
and compare them with experimental values and the results of other first-principles calculations. In 
the liquid-state simulations, we address particular effort to the study of finite-size effects, Brillouin- 
zone sampling and other sources of technical error. Results for the radial distribution function, the 
diffusion coefficient and the shear viscosity are presented for a wide range of thermodynamic states 
relevant to the Earth's core. Throughout this range, liquid iron is a close-packed simple liquid with a 
diffusion coefficient and viscosity similar to those of typical simple liquids under ambient conditions. 
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I. INTRODUCTION 

The aim of this work is to use first-principles simu- 
lation to determine some of the basic properties of liq- 
uid iron at temperatures and pressures relevant to the 
Earth's core, and to present tests and comparisons which 
establish the reliability and robustness of the techniques 
employed. Many lines of evidence show that the Earth's 
core consists mainly of iron, with a minor fraction of light 
impurities |l]-|| . Since the inner part of the core is solid 
and the outer part is liquid, a firm grip on the properties 
of both solid and liquid iron under the relevant conditions 
is crucial to understanding the behaviour of the core, in- 
cluding its convective dynamics, heat transport, and the 
generation of the Earth's magnetic field. 

The laboratory investigation of iron under core con- 
ditions is exceedingly difficult, because of the high pres- 
sures and temperatures needed. At the boundary be- 
tween the mantle and the core, at a depth of ~ 3000 km, 
the pressure is 135 GPa and the core temperature is 
believed to be ~ 4000 K [Q, while at the boundary 
between inner and outer core the pressure is 330 GPa 
and the temperature is thought to be in the region of 
5000 K. Experiments using the diamond anvil cell |^-||] 
can be performed up to ~ 200 GPa, and x-ray diffrac- 
tion measurements on solid Fe have been made at con- 
ditions approaching the core-mantle boundary [||J^] . Be- 
yond this, shock measurements have given some infor- 
mation about the thermodynamic properties of solid and 
liquid iron ||Io|-|l2||, but the temperature calibration of 
these measurements is difficult, and there are substantial 
disagreements between the results. The shear viscosity 
of the liquid is important for understanding core dynam- 
ics, but has been extremely controversial, with estimates 



from different sources spanning many orders of magni- 
tude |l3|,|l4|. The laboratory measurement of viscosities 
using the diamond anvil cell is becoming feasible up to 
pressures of ~ 10 GPa Jl5|,[l6|], but great technical prob- 
lems still need to be overcome before such measurements 
can be done under core conditions. 

First-principles calculations based on density func- 
tional theory (DFT) ||l7|-[l9|] are becoming increasingly 
important in the study of materials under extreme con- 
ditions, and a substantial effort has already been de- 
voted to iron. It is well established that DFT ac- 
curately reproduces the properties of low-temperature 
body-centred cubic (b.c.c.) iron at ambient pressures, 
including the equilibrium lattice parameter, bulk modu- 
lus and magnetic moment ]20| , |2l"| ] , and phonon frequen- 
cies 22 1 . There has also been much DFT work on dif- 
ferent crystal structures of Fe at high pressures, and ex- 
perimental low-temperature results for the pressure as 
a function of volume p(V) up to p = 300 GPa for the 
hexagonal close-packed (h.c.p.) structure are accurately 
predicted |pfj| , pT|j23| . Further evidence for the accuracy of 
DFT comes from the successful prediction of the transi- 
tion pressure from the b.c.c. to the h.c.p. phase p(i| , p3"| |. 
Building on these successes, there have been extensive 
investigations [§2|,^3| of the relative stability and proper- 
ties of the main candidate crystal structures at pressures 
up to core values. 

Very recently, we have applied 

the pseudopotential/plane-wave formulation of DFT to 
calculate the properties of both solid and liquid iron un- 
der core conditions pl| , |24] -p6[| . The great advantage of 
this DFT formulation is that systems containing tens 
or hundreds of atoms can be treated. Furthermore, the 
forces on the atoms are readily calculated, so that dynam- 



1 



ical simulations can be performed, and liquids in thermal 
equilibrium can be simulated. We have reported a brief 
study of the liquid structure , which showed that it is 
close packed, with a coordination number slightly above 
12. We also reported values for the diffusion coefficient, 
and we used the approximate Stokes-Einstein relation to 
infer values of the viscosity (^4|. Our results indicated 
that the viscosity of l-Fe under core conditions is roughly 
10 times greater than that of typical simple liquid metals 
under ambient conditions. However, our liquid Fe sim- 
ulations were performed on a very small system of 64 
atoms, and we were unable at that stage to investigate 
system-size errors or errors due to other technical factors 
such as fc-point sampling. In addition, we studied only a 
very small number of thermodynamic states. 

Here we report the results of much more extensive first- 
principles simulations covering a wide range of thermo- 
dynamic conditions, and we also describe the thorough 
study of technical issues that we have now made. An im- 
portant advance over our previous work is that most of 
the present calculations are based on the projector aug- 
mented wave (PAW) formulation of DFT [|f|j2|. This 
is an all-electron technique which is closely related both 
to other standard all-electron techniques like the linear 
augmented-plane-wave (LAPW) method , and also to 
the ultrasoft pseudopotential method used in our previ- 
ous work pl| , |24] -|26]| . However, in contrast to other all- 
electron methods, PAW allows one to do dynamical first- 
principles simulations of the kind that are routinely done 
with the pseudopotential approach. We report detailed 
tests on different crystal structures of Fe, which show 
that the PAW technique reproduces very accurately the 
available experimental data, as well as the results of pre- 
vious all-electron calculations. We describe the results of 
our liquid-state PAW simulations on a range of different 
system sizes, which allow us to give a quantitative assess- 
ment of size errors and other technical factors. A further 
advance over our earlier work is that we now calculate 
the viscosity directly from the Green-Kubo relation in- 
volving the stress autocorrelation function |50|, rather 
than estimating it from the diffusion coefficient. 

Details of our DFT techniques are given in the next 
Section, and the results of our tests on crystalline iron 
are reported in Sec. 3. Our detailed tests on the reliabil- 
ity of our techniques for studying the liquid are presented 
in Sec. 4, where we also report results for the radial dis- 
tribution function, diffusion coefficient and viscosity of 
the liquid over a range of conditions. Discussion and 
conclusions are given in Sec. 5. 



II. METHODS 

Density-functional theory is a general and extremely 
widely used set of techniques for treating the energetics 
of condensed matter, and it has frequently been reviewed 



(see e.g. Refs. ]18| , |l9| , pl[]32| ). Originally formulated as 
a method for calculating the energy of the electronic 
ground state of collections of atoms for given nuclear po- 
sitions it was later generalized to calculate the 
electronic free energy when the electrons are in thermal 
equilibrium at some finite temperature |35| (again, for 
given nuclear positions). Since DFT also yields the forces 
on the nuclei via the Hcllmann-Feynman theorem, it is 
also possible to perform first-principles molecular dynam- 
ics |l7j , in which the nuclear positions evolve according to 
classical mechanics under the action of the forces, while 
the electronic subsystem follows adiabatically, either in 
the ground state or in the state of instantaneous thermal 
equilibrium. 

Implementations of DFT can be divided into all- 
electron (AE) methods, in which both core and va- 
lence electrons are treated explicitly, and pseudopoten- 
tial methods, in which only valence electrons are treated 
explicitly, their interactions with the cores being rep- 
resented by a pseudopotential. AE methods, such as 
the well-known FLAPW (full-potential linear augmented 
plane-wave) and FP-LMTO (full-potential linear muffin- 
tin orbitals) schemes, are traditionally regarded as more 
'rigorous', but at present can be applied only to rather 
small numbers of atoms, and it is difficult to use them for 
dynamical simulations. The pseudopotential approach, 
by contrast, is routinely used for dynamical simulations 
on systems of tens or hundreds of atoms. First-row and 
transition- metal elements used to be far more demanding 
with the pseudopotential approach, but the introduction 
of so-called 'ultrasoft' pseudopotentials (USPP) |36) has 
largely overcome these problems. 

In the last few years, a new method, known as PAW 
(projector augmented wave) has been developed p7|, 
which effectively bridges the divide between AE and 
pseudopotential methods. It is an all-electron method, 
in the sense that it works with the true Kohn-Sham or- 
bitals, rather than orbitals that are 'pseudized' in the 
core regions, and it has the same level of rigor as AE 
methods such as FLAPW. At the same time, it is very 
closely related to the USPP technique, and reduces to 
this if certain well-defined approximations are made, as 
shown in the recent analysis by Kresse and Joubert [^8) . 
The nuclear forces can be calculated in the same way 
as in pseudopotential methods, so that dynamical sim- 
ulations can be performed at essentially the same cost 
as with the USPP method. An extensive set of com- 
parisons between the PAW, USPP and AE approaches 
has recently been presented [^8| for a variety of small 
molecules and bulk crystals, including Fe, Co and Ni, 
and it is was shown that the three approaches give virtu- 
ally identical results in most cases. The only significant 
exception is ferromagnetic Fe, where fully converged re- 
sults are easier to obtain with PAW than with USPP. The 
present work makes direct use of the work of Kresse and 
Joubert, and was performed with the PAW facility devel- 
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oped in the VASP code [g7|3J]. The details of the core 
radii, augmentation-charge cut-offs etc. in the present 
PAW work are exactly as in the Fe calculations reported 
in Ref. §|. 

It is known that under Earth's core pressures it is not 
accurate enough to neglect the response of the 3s and 
3p electrons to the high compression. The 3p response is 
most significant, and in order to achieve good accuracy 
we treat both 3p and 3d as valence electrons. However, 
this makes it very costly to do long dynamical simula- 
tions on large systems, and we showed earlier |^l| that 
it is an accurate approximation to treat the 3p orbitals 
as rigid provided the effect of their response is included 
by an effective pair potential. (We stress that non-linear 
core corrections are included everywhere in the present 
work.) The procedure used to construct this effective 
potential is as follows (see also Rcf. We calculate 

the total energies of an perfect h.c.p. Fe crystal in two 
ways: in the first we include the 3p orbitals in valence, 
and in the second we freeze them in the core. The effec- 
tive pair potential is then constructed so as to reproduce 
the difference of the energies between the two sets of cal- 
culations. The accuracy of this procedure for treating 3p 
response is demonstrated in Sec. 3 in our comparison of 
the phonon frequencies of the f.c.c. phase at high pres- 
sure. In the following, we refer to calculations in which 
3p and 3s orbitals and below arc in the core as 'Ar-core' 
calculations, and those in which 3p orbitals are included 
in the valence set but 3s orbitals and below are in the 
core as 'Ne3s 2 -core' calculations. 

In the temperature range T > 3000 K of interest 
here, thermal excitation of electrons is important, and 
throughout this work we use the finite-temperature for- 
mulation of DFT |J5|. This means that each Kohn- 
Sham orbital has an occupation number given by the 
usual Fermi-Dirac formula. This thermal excitation has 
an important effect both on the pressure in the system 
and on the nuclear forces, as will be reported in detail 
elsewhere. However, it should be noted that this pro- 
cedure ignores the possible temperature dependence of 
the exchange-correlation functional about which little is 
currently known. 

Further technical details are as follows. All the calcu- 
lations presented here are based on the form of general- 
ized gradient approximation (GGA) known as Perdew- 
Wang 1991 ffl,E3]. In the very few spin-polarized cal- 
culations we have performed, the spin interpolation of 
the correlation energy due to Vosko et al. was used [^2) . 
Brillouin-zone sampling was performed using Monkhorst- 
Pack (MP) special points 43 1. Most of the dynamical 
simulations on the liquid were performed using T-point 
sampling only, though we shall report tests using more 
fc-points. The plane wave cut-off was 300 eV in all calcu- 
lations. The extrapolation of the charge density from one 
step to the next was performed using the technique de- 
scribed by Alfe Q . The time-step used in the dynamical 



simulations was 1 fs. 



III. SOLID IRON 

In order to demonstrate the robustness and reliabil- 
ity of the present PAW methods, we have calculated a 
range of solid-state properties of iron and compared them 
with experimental values, and with the results of other 
DFT implementations, including FLAPW, FP-LMTO 
and ultra-soft pseudopotentials. We also probe the ef- 
fect of treating the valence-core split in different ways. 
We have taken pains to ensure that every result given by 
the present calculations is fully converged with respect to 
fc-point sampling, plane-wave cut-off, and all other tech- 
nical factors. 

In Table ffl, we report values of the equilibrium volume 
Vo for the b.c.c. and h.c.p. structures, the bulk modulus 
Kq and its pressure derivative dKo/dp, and the magnetic 
moment /i per atom of the (ferromagnetic) b.c.c. phase. 
In the calculation of Kq and K' for b.c.c, the moment 
Pl is, of course, free to change with volume. In the h.c.p. 
case, the equilibrium value of the c/a ratio is determined 
by minimizing the total energy at fixed volume. Since in 
principle c/a can depend on volume, this must be taken 
into account in the calculations of Kq and K' . In the re- 
gion of zero pressure we find the value c/a — 1.58, which 
agrees exactly with the calculated value of Ref. (2^] , and 
is in reasonable agreement with the experimental values 
in the range 1.58 — 1.61 reported in Ref. Q.We note 
the excellent stability of the results with respect to DFT 
implementation, and the generally very good agreement 
with experimental values. We expect the most accurate 
variant of our PAW calculations to be the one based on 
the Ne3s 2 core, and Vq values calculated in this way agree 
with FLAPW values to within 1 % or better. The agree- 
ment between PAW[Ne]3s 2 and FLAPW results for Kq 
and Kq is also excellent for h.c.p. . The magnetic mo- 
ment depends only weakly on the DFT method. The 
agreement of our PAW results with experimental values 
is very good for b.c.c, but less good for h.c.p.. However, 
the experimental values for Vq and Kq in the h.c.p. phase 
do not come from direct measurements at zero pressure 
(since the h.c.p. crystal is unstable at this pressure), 
but from an extrapolation from measurements at higher 
pressures. The reliability of the "experimental" values is 
therefore not clear. 

Our calculated phonon dispersion curves for the zero- 
pressure b.c.c. phase are compared with experimental 
frequencies in Fig. ^. The calculations were done using 
our implementation of the small-displacement method 
described in Ref. pq . Small systematic differences can 
be seen for transverse modes at some wavevectors on the 
r — N and T — H lines, but even at worst these do not 
exceed ca. 10 %. 
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Fig. |2| compares calculated values for p as a function of 
V of the h.c.p. crystal with experimental measurements 
due to Mao et al. [Q. On the scale of the Figure, the 
PAW results obtained with the Ne3s 2 core are indistin- 
guishable from those given by the Ar core plus pair po- 
tential, and are also virtually identical to the FP-LMTO 
results of Soderlind et al. p3[ . There is a detectable dif- 
ference from the FLAPW results of Stixrude et al. [ p0[ , 
but this is smaller than the difference from the experi- 
mental values, which itself is very small in the pressure 
range 100 — 300 GPa of main interest in this paper. 

A further test concerns the low-temperature coexis- 
tence pressure of the b.c.c. and h.c.p. phases, which 
experimentally is in the range 10 — 15 GPa |47[| . The 
earlier FLAPW and FP-LMTO calculations both gave 
transition pressures of ca. 11 GPa. The present PAW 
calculations give 10, 12 and 13 GPa for the PAW[Ne]3s 2 , 
PAW[Ar] and PAW[Ar] + effective pair potential respec- 
tively, in reasonable agreement with other values. 

Finally, we report our comparison of phonon frequen- 
cies of the high-pressure non-magnetic f.c.c. phase with 
the FP-LMTO calculations of Soderlind et al. @. Fre- 
quencies were computed at the three volumes and four 
wavevectors studied in the FP-LMTO work, and their 
values are compared in Table |n[ We note that this is not 
a very direct comparison, since the FP-LMTO calcula- 
tions were based on the LDA rather than the GGA used 
here. We have shown elsewhere [M] that the frequencies 
of f.c.c. Fe increase slightly when LDA is replaced by 
GGA, and this is consistent with the difference shown 
in the Table. However, the differences are no more than 
a few percent. We also note that the PAW frequencies 
obtained with the Ar core plus pair potential are almost 
identical to those given by the Ne3s 2 core, the effect of 
the pair potential is to increase the frequencies by 5-10%. 
This is useful evidence that the procedure of using the 
Ar core plus pair potential gives a good account of vibra- 
tional as well as equilibrium properties. 

Our overall conclusion from these tests is that our 
present PAW calculations appear to reproduce very ac- 
curately the properties of real crystalline Fe, and the dif- 
ferences from pseudopotential and standard all-electron 
implementations of DFT are extremely small. This pro- 
vides a firm foundation for our PAW simulations of liquid 
Fe presented in the next Section. 

IV. LIQUID IRON 

A. A representative ab initio simulation 

For orientation purposes, we present first the results 
obtained from a direct ab initio simulation of liquid 
iron at the temperature T = 4300 K and the density 
p = 10700 kg m~ 3 , which corresponds roughly to condi- 
tions at the core-mantle boundary. This simulation was 



performed on a system of 67 atoms, and used T-point 
sampling. The system was initiated and thoroughly equi- 
librated in a way that will be described below. The dura- 
tion of the simulation after equilibration was 15 ps. The 
pressure was calculated in the course of the simulation 
and its value was 132 GPa, which should be compared 
with the value 135 GPa at the core-mantle boundary. 

The radial distribution function g(r) from this simu- 
lation is displayed in Fig. ||. It shows the form typical 
of real and theoretical simple liquids such as liquid alu- 
minum or the Lennard- Jones liquid. It is instructive to 
calculate the average coordination number N c character- 
ising the number of nearest neighbours surrounding each 
atom. This is defined as: 

N c = A"Kn I drr 2 g(r) , (1) 
J o 

where n is the atomic number density and r c is the po- 
sition of the first minimum in g(r). We find the value 
N c = 13.2. This is actually slightly greater than the co- 
ordination number of 12 exhibited by the close-packed 
h.c.p. and f.c.c. structures, and it is clear that the liq- 
uid structure displays a very dense form of packing. We 
note that by no means all monatomic liquids show this 
high coordination. Directionally bonded liquids such as 
/-Se and /-Si have very much lower coordination numbers 
(N c ~ 2 and 6 in these two cases). The high value of 
N c we find for high-pressure /-Fe is evidence of strong re- 
pulsive forces between the atoms and a complete absence 
of directional bonding. This is, of course, not surprising 
given the very high degree of compression. 

We now pass to the self-diffusion coefficient D, which 
we obtain from the asymptotic slope of the time- 
dependent mean-square displacement (MSD) (| Ti(t + 
to)-ri(t ) | 2 }: 

P(t) = (| n{t + 1 ) - n(t ) | 2 ) -> 6D | t | (2) 

in the long-time limit | t |— » oo. The quantity P(t)/6t 
from our simulation is reported in Fig. |]. We also report 
in the Figure the statistical errors on the MSD, which we 
shall need to refer to later; these errors were estimated 
by calculating the MSD separately for every atom in the 
system and examining the scatter of the results. As gen- 
erally happens in simple liquids, the asymptotic linear 
region in the MSD is attained extremely quickly, after a 
transient period of only ~ 0.3 ps. The asymptotic slope 
gives the value D = (5.2 ±0.2) x 10~ 9 m 2 s _1 . The value 
of D is quite typical of simple liquids under ambient con- 
ditions. For example, the diffusion coefficient of liquid 
Al near its triple point is ~ 8 x 10~ 9 m 2 s _1 . This 
indicates that the dynamics of the atoms in /-Fe under 
the present thermodynamic conditions is not drastically 
affected by the high degree of compression. 

The final quantity we examine is the shear viscosity 
rj. In our earlier work p4j |, we relied on the Stokes- 
Einstein relation to give a rough estimate of the viscosity, 
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but we subsequently showed that this indirect method is 
completely unnecessary, since the direct, rigorous calcu- 
lation of rj using the Green-Kubo relation is perfectly 
feasible |p0| . According to this relation, 77 is given by: 

where (P X y(t)P X y(0)) is the stress autocorrelation func- 
tion (SACF), i.e. the autocorrelation function of the off- 
diagonal component P xy of the stress tensor at times 
separated by t, and V is the volume of the system. 
Techniques for calculating the SACF are presented in 
Rcf. j30|, where it is noted that statistical accuracy is 
improved by taking the average cj>(t) of the five indepen- 
dent correlation functions of the traceless stress tensor 

Pxy i Pyz? Pzx: "n\Pxx ^yy) &ud "n(Pyy Pzz}- 

The average SACF (f>(t) and its time-integral 
J d£ 4>(t') for the present thermodynamic state are re- 
ported in Fig. H with error bars for the time integral. 
We obtain from this the value 77 = 8.5 ± 1 mPa s. It is 
interesting to compare this with the value that would be 
obtained from the Stokes-Einstein relation: 

Dr] = k B T/2Tra , (4) 

where a is an effective atomic diameter. In our previous 
work, we took the effective atomic diameter a to be the 
nearest-neighbour distance in the close-packed solid hav- 
ing the same density, which in the present case is 2.15 A. 
Using the value D = 5.2 x 10 -9 m 2 s _1 reported above, 
we obtain the estimate n = 8.5 mPa s, which happens to 
be the same as the Green-Kubo value. 

Useful though they are, the results we have just pre- 
sented are limited in two ways. First, there are technical 
limitations. Our results have been obtained from a rather 
small simulated system of only 67 atoms, and we must 
clearly try to show that they are not seriously influenced 
by size effects. Similarly, we need to know that the re- 
sults are not affected by other technical factors, such as 
the use of T-point sampling, the choice of the PAW tech- 
nique rather than some other ab initio technique, or the 
split that we have chosen between valence and core states. 
The second limitation is that we have studied only one 
thermodynamic state. We want to know how the prop- 
erties of l-Fe vary over the range of conditions relevant 
to the Earth's core. In overcoming both these kinds of 
limitation, the notion of a reference system will be very 
important, and we discuss this next. 

B. The reference system 

Since ab initio simulation is very costly, it is difficult 
to study size effects by directly simulating large systems. 
This is particularly difficult for dynamical properties like 
the diffusion coefficient and the viscosity, since rather 



long simulations are needed. It will be therefore very 
helpful to have available a simple model or reference sys- 
tem which closely mimics the behaviour of the full ab 
initio system. We have in mind a model in which the 
atoms interact through a simple pairwise potential (f>(r), 
so that its total potential energy Uq is given by: 

If the reference system is to behave like the ab initio sys- 
tem, then its total energy Uo must closely resemble the 
ab initio total energy U. This means that the difference 
AU = U — Uo must be small. In fact, the requirement is 
slightly weaker than this. It makes no difference to the 
thermal-equilibrium structure or dynamics of the liquid if 
we add a constant to the total energy. So the requirement 
is really that the fluctuations of AU should be small. We 
do not need to demand that the strength of these fluc- 
tuations be small over the whole of configuration space. 
It suffices that they are small over the region explored 
by the atoms. Regions of configuration space rarely vis- 
ited by the system should be given a low weighting in 
assessing the fluctuations, and the natural way to do this 
is to weight configurations with the Boltzmann factor 
exp[— (3U(ri, . . . rjv)]- The requirement is therefore that 
the thermal average of the mean-square fluctuations of 
AU defined by ((AU — (AU)) 2 } be as small as possi- 
ble. Here, the thermal average ( ■ ) can be evaluated as 
a time average over configurations generated in the ab 
initio simulation. It will be noted that this requirement 
depends on thermodynamic state, so that the best model 
system may be different for different states. 

We have every reason to expect that our ab initio liquid 
can be well modelled by a simple reference system, since 
we have seen that its properties closely resemble those 
of typical simple liquids. Initially, we attempted to use 
the Lennard- Jones model, in which cf>(r) = 4e((a/r) 12 — 
(r/cr) 6 ), where a and e characterise the atomic diame- 
ter and the depth of the attractive potential well respec- 
tively. We found that it is possible to choose a and e 
so that the fluctuations SAU = AU — (AU) are reason- 
ably small: the minimum rms value that we achieved 
was [((SAU) 2 )^] 1 ' 2 = 0.18 eV (N is the total number 
of atoms in the system). However, the rdf of the result- 
ing LJ liquid differs appreciably from that of the ab initio 
system. Further experimentation showed that a far bet- 
ter reference system is provided by a pure inverse-power 
potential <f)(r) — B/r a . After adjusting the parameters 
B and a so as to minimise the strength of the fluctu- 
ations SAU, we achieved a much improved rms value 
[((SAU) 2 )^] 1 / 2 equal to 0.08 eV. Bearing in mind that 
ksT = 0.37 eV at the temperature of interest, this repre- 
sents a very satisfactory fit to the ab initio system. The 
resulting value of the exponent a is 5.86, and B is such 
that for r = 2 A the potential <j>(r) is 1.95 eV. 
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The inverse-power reference system is so close to the 
ab initio system that its rdf go (r) is almost indistinguish- 
able from the rdf of the ab initio system. (Naturally, in 
making this comparison, we treat the two systems at ex- 
actly the same thermodynamic state and using the same 
67-atom repeating cell.) To show this, we plot the differ- 
ence Ag(r) = g[f) — go(r) in Fig. ||. The first peak has 
almost exactly the same height in the two cases, and the 
first minimum is also virtually identical; the difference of 
the two rdfs consists mainly of a slight shift of the first 
peak of the reference-system to a smaller radius. 

Since the reference system is so good, we should expect 
it to reproduce also the dynamics of the ab initio system. 
To test this, we compare in Fig. ^ the msd of the ref- 
erence system with the ab initio results, using again the 
67-atom repeating system. The two curves agree very 
closely for short times, but there is a significant difference 
for t > 0.1 ps, which is well outside the statistical errors. 
The diffusion coefficients obtained for PAW and the ref- 
erence model are 5.2 x f 0~ 9 and 6.f x 10 -9 m 2 s _1 respec- 
tively. Given that the peak position of g(r) is at a slightly 
smaller radius for the reference model, it is perhaps sur- 
prising that the model makes the atoms more mobile, 
but one should remember that the rdf describes only 2- 
body correlations, and higher correlations may well be 
important for diffusion. 

Finally, we make the same comparisons for the stress- 
stress correlation function and its time integral (see 
Fig. ||). The asymptotic value of the viscosity integral 
is significantly smaller for the reference model - as one 
would expect from the higher diffusion coefficient in this 
case. The ab initio and reference viscosities are 8.5 and 
7.0 mPa s respectively, so that they differ by ca. 20 %, 
as would be expected from the Stokes-Einstein relation. 

Our overall conclusion from these comparisons is that 
the simple inverse-power reference system with appropri- 
ately chosen parameters reproduces the structure of the 
ab initio liquid very well and its dynamics reasonably 
well. It will be asked whether the same reference system 
also works for other thermodynamic states, or whether 
the parameters B and a have to be refitted for every 



A duration of 1 ps is enough to give excellent statistical 
accuracy for the rdf, so our examination of size effects on 
g(r) was done using the ab initio simulations directly. 
We find that the dependence of g(r) on system size is so 
small that it cannot easily be seen on simple plots of g(r). 
To show these effects, we therefore plot the differences 
5A r ( r ) — 9N'(r) between the rdfs for different numbers 
of atoms. These differences are reported in Fig. || for 
(A, A') equal to (127,67) and (241,127). The difference 
between g^-j (r) and 5127 (r) is small and comes almost 
entirely in the region of the first peak; the same is true 
of the difference between (7127 (r) and 3241 (r). 

We have also looked for size effects by studying how 
well the inverse-power reference system constructed by 
fitting to the 67-atom system reproduces the total energy 
of the larger systems. We find that the rms strength of 
the fluctuations [((SAU) 2 ) /A] 1 / 2 observed in the 1 ps ab 
initio simulations for the larger systems is essentially the 
same as what we found for the 67-atom system. This 
confirms that the reference system mimics the large ab 
initio systems as well as it mimics the small one. 

A duration of 1 ps is too short to study size effects on 
D and n directly from the ab initio simulations: com- 
parisons between results for different system sizes would 
be vitiated by statistical noise. However, since the refer- 
ence model appears to work equally well for all system 
sizes, we can legitimately use this model to study size 
effects. We have therefore calculated D and 77 from sim- 
ulations of the reference model, using simulations lasting 
for 500 ps after equilibration. The results are reported 
in Table III. We note that D increases slightly as we 
go to larger systems. The effect of system size on D 
in hard-sphere systems was studied by Erpenbeck and 



state. This question will be answered in Sec. IV E 



Wood 1 49 1, who showed that the calculation of D with 
64 atoms would underestimate its value by ca. 20%. This 
is in the same direction as the effect we are seeing, but 
in our case the error appear to be less than 10%. Ac- 
cording to the Stoke-Einstein relation, we would expect 
a decrease of rj with increasing system size, but this is 
not clear from our results. 



D. Other technical factors 



C. Size effects 

As an aid to addressing size effects, we have performed 
a number of ab initio simulations at the same thermody- 
namic state as before (T = 4300 K, p = 10700 kgm~ 3 ), 
using cell sizes ranging from 89 to 241 atoms. The prepa- 
ration and equilibration of these systems were done using 
the inverse-power reference system. Since the latter so 
closely mimics the ab initio system for the 67-atom cell, 
it should provide a well equilibrated starting point for ab 
initio simulation of larger systems. The duration of all 
the ab initio simulations after equilibration was 1 ps. 



We made tests to assess the influence of Brillouin- 
zone sampling. Since ab initio simulations using many 
fc-points are very demanding, we have adopted an in- 
direct approach. Instead of performing direct ab initio 
simulations with many fc-points, we have drawn sets of 
atomic positions {r^} from the T-point simulation, and 
calculated the total energy Up with 4 MP fc-points for 
these configurations. We then compute the rms strength 
of the fluctuations [((5AU P ) 2 ) /A] 1 / 2 , where SAU P de- 
notes the fluctuation Up — U\ — (Up — U%), with U\ the 
T-point energy. This rms fluctuation strength is 0.02 eV, 
which is much smaller then the rms fluctuation between 
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the ab-initio (T-poiiit) system and the reference system. 
This implies that the expected error in D and r\ would 
be barely noticeable, certainly smaller than the statistical 
errors on these quantities. 

In order to test the reliability of using the Ar core plus 
effective pair potential (see Sec. [□]), we did a simulation 
of the liquid at the same thermodynamic state as before, 
but with the Ne3s 2 core. The duration of this simulation 
was 4.5 ps. The rdf for this simulation was almost identi- 
cal to that obtained with the Ar core, the difference being 
somewhat smaller than the size-effect difference between 
g(r) for 67 and 127 atoms discussed above. Finally, we 
have compared g(r) obtained in the present PAW simu- 
lations with the one given by our earlier simulations j2j| ] 
at the same thermodynamic state, which were based on 
ultrasoft pseudopotentials. Here again, the difference is 
no bigger than the size effect between 67 and 127 atoms. 

Our conclusion from these tests is that our simula- 
tions are completely robust with respect to size effects, k- 
point sampling, the split between valence and core states, 
and the use of PAW rather than the pseudopotential ap- 
proach. 



E. Dependence of liquid properties on 
thermodynamic state 



In addition to the simulations reported in Sec. IV A 



we 

have also performed PAW simulations at the 15 thermo- 
dynamic states listed in Table [j^, all these simulations 
being done on the 67-atom system. With these simula- 
tions, we cover the temperature range 3000 — 8000 K and 
the pressure range 60 — 390 GPa, so that we more than 
cover the range of interest for the Earth's liquid core. The 
Table reports a comparison of the pressures calculated in 
the simulations with the pressures deduced by Anderson 
and Ahrens from a conflation of experimental data. 
Our first-principles pressures reproduce the experimental 
values to within 2 — 3 % at low densities, but they are 
systematically too high by ca. 7 % at high densities. It is 
not clear yet whether the high-density discrepancies indi- 
cate a real deficiency in the ab initio calculations rather 
than problems in the interpretation of the experimental 
data. We are currently using free-energy calculations to 
study the thermodynamics of the liquid in more detail, 
and we hope that this will shed light on the question. 

Rather than reporting detailed results for g(r), D and 
•q at each thermodynamic state, we can exploit the prop- 
erties of the reference system to present the results in 
a compact form. We have taken the reference system 
fitted to the ab initio simulations at T = 4300 K and 
p = 10700 kg m~ 3 , and examined the fluctuations 8AU 
in the ab initio simulations done at the other thermody- 
namic states. We find that the strength of the fluctu- 
ations is not significantly greater for these other states 
than it was for the state for which the model was fitted. 



This strongly indicates the same reference model, with 
exactly the same parameters, mimics the first-principles 
system (almost) equally well at all thermodynamic states. 

The invariance of the reference system implies a truly 
remarkable simplicity in the variation of the liquid prop- 
erties with thermodynamic state. It has long been recog- 
nised that the properties of an an inverse-power system in 
thermal equilibrium depend non-trivially only on a single 
thermodynamic variable, rather than on temperature T 
and number density n independently [pl|. This variable 



C = Bn a/3 /k B T 



(6) 



Any other dimensionless combination of the parameters 
B, T and n is a function only of C- This means, for 
example, that the rdf g{r), since it is dimensionless, can 
depend only on £, provided it is regarded as a function of 
the dimensionless distance £ = n x ^r. To formulate this 
statement precisely, let us define the 'reduced' rdf g(£) 
as: 



g(n^ 3 r) = g(r). 



(7) 



Then g(£) must be a function only of £■ Similarly, the 
reduced diffusion coefficient D and viscosity fj, defined 
by: 

D = {m/k B T) 1/2 n 1/3 D , 77 = {mk B T)- 1/2 n- 2/3 r] 

(8) 

can depend only on £. 

To show the independence of the rdf, we report in 
Fig. [7] the rdf g(r) for the five temperatures 4300, 
5000, 6000, 7000 and 8000 K at the same density p = 
10700 kg m~ 3 ; the corresponding values of £ are 4.51, 
3.88, 3.23, 2.77 and 2.43. The effect of varying temper- 
ature (or C) is clearly not dramatic, and consists of the 
expected weakening and broadening of the structure with 
increasing T. Since the variation is so regular, we can ob- 
tain for any intermediate value of C by simple linear 
interpolation between the curves. With this information, 
we can then test explicitly our expectation that g(£) de- 
pends only on £. To do this, we go to states at other 
densities p and compare the reduced rdf with the <7~(£) 
for the same £ but for the density p — 10700 kg m~ 3 . 
We find that the agreement is excellent; indeed, the dif- 
ferences between <?(£) at the same £ but different density 
are even smaller than those between g(r) for ab initio and 
reference model discussed in Sec. IV A. This shows that 



the rdf for any thermodynamic state in the range we have 
studied can be obtained directly from the results shown 
in Fig. |. 

Since most of the simulations reported in Table ^ are 
rather short - typically no more than 4 ps - the statisti- 
cal accuracy on D and 77 is not great. Numerical results, 
together with error estimates obtained as in Sec. IV A, 



7 



are reported in Table [v|. If the reduced quantities D and 
fj defined in eqn (||) depend only on £, then a plot of D 
against £ should consist of a single curve, and similarly 
for fj. This hypothesis is tested in Fig. ||, and appears to 
be satisfied within the (significant) errors. This means 
that the diffusion coefficient and viscosity at any ther- 
modynamic state in the range we have studied can be 
deduced from the data in Fig. ||. 



V. DISCUSSION AND CONCLUSIONS 

The ab initio simulations we have presented demon- 
strate that under Earth's core conditions l-Fe is typical 
simple liquid. In common with other simple liquids like 
l-Ai and l-A\, it has a close-packed structure, the coordi- 
nation number in the present case being ca. 13. For the 
entire pressure-temperature domain of interest for the 
Earth's outer core, the diffusion coefficient D and viscos- 
ity r] are comparable with those of typical simple liquids, 
D being ca. 5 x 1CP 9 m 2 s _1 and rj being in the range 
8 — 15 mPa s, depending on the detailed thermodynamic 
state. 

We have shown that the structural properties of l-Fe 
are reproduced very accurately and the dynamical prop- 
erties fairly accurately by an inverse-power model, and 
consequently that they display a remarkable scaling prop- 
erty. Instead of depending on T and n separately, their 
only non-trivial dependence is on the combined thermo- 
dynamic variable £ discussed in Sec. [VE. Since the 



Earth's outer core is in a state of convection, the tem- 
perature and density will lie on adiabats. It is straight- 
forward to show that the variation of C along adiabats 
will be rather weak. For example, if we take the data 
for high-pressure liquid iron compiled by Anderson and 
Ahrens |o|, then the adiabat for T = 6000 K at the ICB 
pressure 330 GPa has T = 4300 K at the CMB pressure 
of 135 GPa. Taking the densities at these two points from 
the same source, we find values of C equal to 4.85 and 4.51 
at the ICB and CMB respectively. Fig. || shows that this 
small variation of £ leaves the liquid structure virtually 
unchanged, apart from a trivial length scaling. For all 
practical purposes, then, it can be assumed that varia- 
tion of thermodynamic conditions across the range found 
in the core has almost no effect on the liquid structure. 
For the same reasons, the reduced diffusion coefficient D 
and viscosity fj also show little variation. From the results 
reported in Sec. IV E, it is readily shown that the true 
(unreduced) diffusion coefficient D is 5±0.5 x 10~ 9 m 2 s _1 
without significant variation as one goes from ICB to 
CMB pressures along the above-mentioned adiabat, and 
r\ goes from 15 ± 5 mPa s to 8.5 ± 1 mPa s along this 
adiabat. 

We have made strenuous efforts to demonstrate the ro- 
bustness and reliability of the calculations. Our results 



on the solid show that the predictions from our pseudopo- 
tential and PAW ab initio techniques agree very closely 
both with experiment and with the results of other all- 
electron calculations. This is true both at ambient pres- 
sures and at Earth's core pressures. For the liquid, we 
have shown that our results are completely robust with 
respect to all the main technical factors: size of simulated 
system, fc-point sampling, choice of ab initio method and 
split between core and valence states. At present, ab ini- 
tio dynamical simulations of the type presented here are 
only practicable with the PAW or pseudopotential tech- 
niques. However, our comparisons leave little doubt that 
if such simulations were feasible with other DFT tech- 
niques, such as FLAPW, virtually identical results would 
be obtained. In addition, we have been able to compare 
the calculated pressure in the liquid with values deduced 
from experiments, and again the good agreement sup- 
ports the validity of our techniques. 

Finally, we note the implications of this work for the 
controversy about the viscosity of the Earth's core. In 
our earlier work on this problem, we presented ab initio 
simulations on a 64-atom system and used results for the 
diffusion coefficient together with the Stokes-Einstein re- 
lation to argue that the viscosity of pure l-Fe under core 
conditions is ca. 13 mPa s. Although we regarded this 
prediction as sound, it was certainly open to the objec- 
tion that the simulated system was too small, and that 
the viscosity was inferred indirectly. We believe that the 
present work has completely overcome both these objec- 
tions, as well as giving results for the viscosity over a 
wide range of conditions. We fully confirm our earlier 
conclusion that the viscosity of pure liquid iron under 
Earth's core conditions is at the lower end of the range 
of proposed values. 

In conclusion, l-Fe under Earth's core conditions is a 
simple liquid with a close-packed structure and a coordi- 
nation number of 13, and a viscosity going from 15 mPa s 
to 8 mPa s as conditions go from the the inner-core 
boundary to the core-mantle boundary. These results 
are completely robust with respect to system-size effects, 
choice of ab initio method, and other technical factors. 
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n.c.p. PAW [NeJJs 


10.25 
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PAW [ArJ 


10.19 


289 
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PAW [Ar] + corr 


10.26 


288 


4.5 




USPP [Ne]3s 2 


10.35 
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4.4 




USPP [Ar] 


10.32 


281 


4.2 




USPP [Ar] + corr 


10.38 


296 


4.2 




FLAPW 


10.20 


291 


4.4 




expt. 


11.2 


208 







TABLE I. Properties of the b.c.c. and h.c.p. crystal struc- 
tures of Fe at zero pressure obtained in the present PAW and 
USPP calculations compared with FLAPW and experimen- 
tal values. Properties reported are: the equilibrium volume 
per atom Vb, the bulk modulus Kq, the pressure derivative 
Kg = dKo /dp and the magnetic moment fi per atom (units of 
Bohr magneton /is). FLAPW results are those of Ref. [ pc| . 
Experimental results all come from Ref. Q, except for p,, 
which comes from Ref. 15311. 
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L[100] 


T[100] 
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T[lll 
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FP-LMTO 


19.6 


13.6 
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9.56 




PAW [Ne]3s 2 


20.4 
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22.9 


9.61 




PAW [Ar]±corr 


20.2 


13.7 


22.9 


9.42 


7.55 


FP-LMTO 


13.4 


10.1 


15.9 


7.19 




PAW [Ne] 


14.2 


10.5 


16.4 


7.16 




PAW [Ar]±corr 


14.0 


10.3 


16.4 


7.15 


9.70 


FP-LMTO 


6.33 


6.37 


8.93 


4.89 




PAW [Ne]3s 2 


7.23 


6.43 


9.52 


4.65 




PAW Ar]±corr 


7.35 


6.49 


9.64 


4.74 



TABLE II. Zone boundary phonon frequencies (THz units) 
of f.c.c. Fe from the present PAW calculations compared with 
the FP-LMTO results of Ref. Q. Results are shown for 
three values of the volume per atom Vb. PAW results are 
given for both the Ne3s 2 core and the Ar core plus effective 
pair-potential correction (see text). 
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67 


127 


241 


499 


L ) PAW [Ne]3s 2 


4.8 ±0.2 








PAW [Ar] + corr 


5.2 ± 0.2 








Ref. model 


6.1 


6.2 


6.4 


6.5 


PAW [Ne]3s 2 


5.5 ± 2 








PAW [Ar] + corr 


8.5 ±1 








Ref. model 


7.0 


7.5 


7.0 


6.5 



TABLE III. Comparison of ab-initio and reference-model 
values of the diffusion coefficient D and the viscosity r\ for 
simulated /-Fe systems containing 67, 127, 241 and 499 atoms 
at T = 4300 K and p = 10700 kg m -3 . PAW results are 
given for the Ne3s 2 core and the Ar core with pair-potential 
correction (see text). 
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T(K) 
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7000 




161 


181 


264 


375 










(250) 


(350) 


8000 




172 


191 


275 


390 












(360) 



TABLE IV. Pressure (GPa units) calculated in the full set 
of 16 ab-initio simulations of i-Fe. Experimental values (in 
parenthesis) are from Ref. fcol. 







p(k 


rm 3 ) 




T(K) 


9540 


10700 


11010 


12130 


D (10" a m A s _i ) 3000 


4.0 ±0.4 








4300 




5.2 ±0.2 






5000 




7.0 ±0.7 






6000 


14 ± 1.4 


10 ± 1 


9 ±0.9 


6 ±0.6 


7000 




13 ±1.3 11 ±1.1 


9 ±0.9 


r] (mPa s) 3000 


6±3 








4300 




8.5 ±1 






5000 




6±3 






6000 


2.5 ±2 


5±2 


7±3 


8±3 


7000 




4.5 ±2 


4±2 


8±3 



13300 



15 ±5 
10 ±3 



TABLE V. The diffusion coefficient D and the viscosity r\ 
from ab-initio simulations of l-Fe at a range of temperatures 
and densities. The error estimates come from statistical un- 
certainty due to the short duration of the simulation. 



10 



n TaaOl 

1 LHH U J 


TaOOl 


rririril 
LHHHJ 


--Vf-\ -- 


//□ 








□ V \ 

- - ..A.A.. 


□\\\ 








/ 





2 - 



1 - 



N 



H 



1 

Ab-initio 
Ref. system 
Difference 



r(A) 



FIG. 1. Phonon dispersion curves of ferromagnetic b.c.c. 
Fe at zero pressure along the [100], [110] and [111] directions. 
Curves show calculated phonon frequencies, open squares are 
experimental data of Ref. 0] . 



FIG. 3. Radial distribution function g(r) of /-Fe at 
T = 4300 K and p = 10700 kg m" 3 . Solid and dashed curves 
are ab-initio and reference-model results, dotted curve is the 
difference of the two. 
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FIG. 2. Pressure as a function of atomic volume of h.c.p. 
Fe. Solid circles are present PAW calculations; dotted and 
chain curves are FLAPW and FP-LMTO results of Ref. Q 
and [23] respectively, diamonds are experimental values of 
Ref. [45 1 
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FIG. 4. Mean-square displacement P(t) divided by 6t for 
/-Fe at T = 4300 K and p = 10700 kg m~ 3 . Solid and 
short-dashed curves show results for ab-initio and reference 
system, with statistical error shown as dots and long dashes. 
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FIG. 5. (a) Stress autocorrelation function (see eqn. 
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J)) normalised to its t = value of liquid Fe at T = 4300 K and 
p = 10700 kg m -3 . Solid and short-dashed curves show ab-initio and reference-system results, (b) Viscosity integral (see text) 
of l-Fe at the same thermodynamic state. Solid and short-dashed curves as in panel (a); dotted and long-dashed curves show 
statistical errors. 
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FIG. 6. Comparison of radial distribution function g(r) 
from ab-initio simulations of l-Fe using different system sizes. 
Solid, long-dash and short-dash curves show g(r) for systems 
of 67, 127 and 241 atoms; dotted and chain curves show differ- 
ences gN(r)-g N > for (N,N') equal to (127,67) and (241,127) 
respectively. 
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FIG. 7. Variation of g(r) with temperature from ab-initio 
simulations of l-Fe at the fixed density p — 10700 kg m~ 3 . 
Results are shown for the five temperatures T = 4300, 5000, 
6000, 7000 and 8000 K. 
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FIG. 8. Reduced diffusion coefficient D and viscosity 7/ as a function of reduced state variable ( (see text) from ab-initio 
simulations of l-Fe at 12 thermodynamic states. Error bars show statistical uncertainty on each point. 
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